Disorder in a quantum spin liquid: flux binding and local moment formation 
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We study the consequences of disorder in the Kitaev honeycomb model, considering both site 
dilution and exchange randomness. We show that a single vacancy binds a flux and induces a local 
moment. This moment is polarised by an applied field h: in the gapless phase, for small h the local 
susceptibility diverges as x{h) ~ ln(l//i); for a pair of nearby vacancies on the same sublattice, 
this even increases to x{h) ~ ^/i.h\[n{l/hyf'^'^). By contrast, weak exchange randomness does not 
qualitatively alter the susceptibility but has its signature in the heat capacity, which in the gapless 
phase is power law in temperature with an exponent dependent on disorder strength. 

PACS numbers: 75.10.Jm, 75.50.Mm, 75.54.Cx 
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Disorder in many particle quantum systems is not only 
unavoidable but also an invaluable probe. Isolated im- 
purities can generate new phenomena, as in the Kondo 
effect [H . They can reveal elusive properties of their host, 
such as the order-parameter symmetry in a superconduc- 
tor And in some materials a residual defect concen- 
tration may dominate low temperature properties 

Quantum magnets provide a setting for such problems 
and some of the most challenging questions appear in 
systems without conventional long range order. Results 
for the antiferromagnetic spin-half Heisenberg chain il- 
lustrate the possibilities that can arise. In the absence of 
disorder its susceptibility remains finite in the low tem- 
perature limit. Site dilution creates free chain ends J4| 
and leads to a Curie contribution to the susceptibility [5| , 
while disorder in exchange strength generates a random 
singlet phase Id], also with a divergent low temperature 
susceptibility [7|. In higher dimensions a single vacancy 
in a gapped quantum paramagnet again produces a local 
moment and a Curie-like response [8|. Experimentally, 
local probes such as nuclear magnetic resonance can sep- 
arate bulk and defect contributions to susceptibility. In 
particular, studies of some geometrically frustrated mag- 
netic materials that are candidate spin liquids find sub- 
stantial defect contributions [1, Q . 

Opportunities for an exact treatment of such problems 
are rare except in one dimension. We show in this pa- 
per, however, that the features that make the pure Ki- 
taev honeycomb model [l3| tractable extend to a system 
with site dilution or exchange randomness. Without dis- 
order this model provides a solvable instance of a spin 
liquid, with both gapped and gapless phases, which sup- 
port fractionaliscd excitations, abelian and non-abelian 
[lO| . It has extensions to a lattice where the exact ground 

to higher spin 



state is a chiral spin liquid 
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to higher dimensions [l3|. In addition, there arc propos- 
als for experimental realisations using cold atoms in an 
optical lattice 14 [ , and as a solid state system [l^ . 

The solvability of the Kitaev model rests on the exis- 
tence of a set of non-dynamical fluxes of an emergent Z2 



gauge field. Within each flux sector the Hamiltonian can 
be reduced to a free fermion problem. Disorder in the 
form of site dilution or exchange randomness is not an 
obstacle to this reduction but has dramatic effects on the 
physical behaviour. 

The simplest form of disorder involves a single va- 
cancy. We find that the vacancy binds a flux. On 
the sites adjacent to the vacancy a local moment forms 
with a signature in the susceptibility. For the pure 
system the ground-state susceptibility is finite [loj . In 
contrast, we show that the vacancy susceptibility is di- 
vergent, varying at weak field h in the gapless phase 
as x(/i) ^ ln(l//i). The moments of different vacan- 
cies interact. Strikingly, two nearby vacancies on the 
same sublattice have a greatly enhanced susceptibility: 
x{h) ~ (/i[ln(l//i)]'^/^)~^. Our results complement re- 
cent work on magnetic impurities in this model [lij . 

We also study the effects of weak randomness in ex- 
change interactions. This type of disorder docs not qual- 
itatively influence the susceptibility, which is left finite. 
It docs however change the form of the heat capacity 
in the gapless phase, since the associated free fermion 
problem is represented at long wavelengths by a massless 
Dirac Hamiltonian with a random vector potential. The 
resulting density of states is power-law in energy [l7j |. 
giving a low-temperature heat capacity that is power- 
law in temperature, with exponents that in both cases 
are continuously dependent on disorder strength. 

The Kitaev honeycomb model is defined as follows. A 
bond on this lattice between nearest neighbour sites j 
and k takes one of three orientations, labelled using the 
variable ajk = x, y or z as shown in Fig. [Ija). A spin 
one-half variable at site j is represented by the Pauli 
matrices a". Exchange of strength Ja^^ acts between 
spin components ajk and the Hamiltonian is 



{Jk) 



(1) 



It is exactly soluble, for example, with a local transfor- 
mation cr" = ib°'c which represents each spin using the 
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four Major ana fermions b^^b^ , b^ and c [l^l- Then 

^ ^-jYl Jo,,^^okCjCk Ujk = ib'^^'b^'" . (2) 

The operators Ujk commute with each other and with 
H. One can therefore fix the values of {iijk) = Ujk = ±1, 
move to a subspace of the fuh Hamihonian and obtain a 
bihnear form in the Cj's. Numbering the sites around a 
plaquette from 1 to 6 [see Fig. [IJa)] the Z2 flux through 
a plaquette is defined to be Wp = M21W23M43M45M65U61. 
Physical properties of the system depend only on these 
fiuxes [10| but note that many choices of the set {ujk} 
give rise to the same flux sector. The ground state sector 
is flux free, e.g. with all ujk = +1 for sites j on a partic- 
ular sublattice [l^l- There is, however, a complication: 
in transforming to Majorana fermions, the Hilbert space 
per spin is doubled and projection is necessary to obtain 
physical states of the system Nevertheless, it can be 
shown [l9| that matrix elements evaluated in a subspace 
with a flxed set of {ujk} are the same as those obtained 
using the projected physical states. 

In a given flux sector Eq. ^ has the form 
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Here, for a lattice of N unit cells, ca and cb are each N- 
component vectors of Majorana fermion operators, from 
sublattices A and B respectively, while M is an x iV 
matrix, with entries UjkJajk- This Majorana Hamilto- 
nian has the same ground state energy as the complex 
fermion Hamiltonian 

"){::) « 

at half filling, and we base our study on the associated 
tight binding model. In the parameter space with all Jq, 
non-negative there are three gapped phases (where one 
exchange dominates: Jz > Jx + Jy or permutations) and 
one gapless phase, around the point Jx=Jy=Jz^J on 
which we focus below. 

To study moment formation due to a vacancy, we intro- 
duce a magnetic field h with Zeeman coupling to spins. 
This contribution to the Hamiltonian renders flux dy- 
namical and so spoils solvability As gaps from 
the ground state to other flux sectors are 0{J), at field 
strength h <Si J one can project onto a given flux sector 
and work perturbatively in h. For the undiluted lattice, 
matrix elements of the Zeeman energy between states 
from the same flux sector are all zero and the leading 
term in a projected perturbation is second order in h 
[lo| . In the presence of a vacancy, however, the individ- 
ual fluxes through the surrounding three plaquettes lose 
physical significance and only their Z2 sum enters the 
definition of a flux sector. Because of this, the Zeeman 




(b) (c) 

FIG. 1; (a) Labelling x,y,z for bond orientations; site num- 
bering 1...6 for plaquette operator Wp. (b) Sites at which 
field components hx, hy, and hz contribute finearly to the pro- 
jected Hamiltonian with a vacancy, (c) New sites of the equiv- 
alent tight binding model, coupled with hopping hx,hy,hz. 



energy acquires non-zero matrix elements within a sector, 
and the first terms in a projected Hamiltonian arc now 
linear in h. These terms arise from the specific field com- 
ponents and sites indicated in Fig. mb): employing the 
site labelling shown there, the projected Zeeman energy 
is Hz = —{hx<Ji -\- hya\ + h^a^). Written using Ma- 
jorana fermions, the contributions to Hi have the form 
hauf = ihabfci. They can be represented in the tight 
binding model, Eq. (|3]), by the addition of new sites, with 
the coupling shown in Fig.[ljc). 

The task now is to calculate the field-dependence of the 
ground state energy of the Kitaev model with a vacancy, 
using the tight binding model of Fig.[TJc). Consider first 
the behaviour in a gapped phase. A consequence of the 
bipartite structure of this tight binding model, in any 
flux sector, is that finite energy eigenstatcs appear only 
as positive and negative energy pairs. There may in ad- 
dition be zero energy states. At h = there are four of 
these, but if all components of h and all Ja are non-zero 
there are only two such states, accompanied by a finite 
energy pair with a separation that is linear in h for small 
/i. As a result, the ground state energy decreases linearly 
with h, reflecting a flnite vacancy moment which varies 
continuously with the exchange parameters at vanishing 
/i in a gapped phase. 

Vacancy properties in the gapless phase arc more sub- 
tle, and cannot be found by considering only a finite num- 
ber of energy levels. We instead proceed by computing 
the Green function G(h) for this tight binding model. It 
is a function of a complex energy parameter z. Defining 
£[z,h) = zTr [G(h) — G(0)], we evaluate the change in 
ground state energy within the flux sector for the Kitaev 
model with a vacancy and field, compared to its zero field 
value, as an anti-clockwise integral around the negative 
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real axis: 

£(h) = { £(z,h)dz . (5) 
27ri J 

As advertised above, the presence of a vacancy changes 
the flux sector in which the ground state lies. However, 
we postpone discussion of this aspect and focus initially 
on the simpler behaviour in the flux free sector, which we 
will make use of for the vacancy pair calculation below. 

In essence, our calculation of £{z,h) depends on the 
fact that G'(h) is related by a finite rank perturba- 
tion to the Green function Go for the undiluted lattice, 
which is known [2^. Moreover, at small h the integral 
in Eq. ^ is dominated by contributions from small z, 
and so results are governed by the behaviour of Gq for 
\z\ ^ J. The standard T-matrix formalism casts G(h) 
in terms of the matrix elements Go(r, z) of Go between 
sites with separation r. These enter in the combination 
g{z) = Go(0, z) -I- Go(l, 2:)^/Go(0, 2), which is the site- 
diagonal element of G(0) at one of the sites (e.g. ri) 
adjacent to the vacancy. Expressions are simplest if only 
one component of h is non-zero, and in this case 

£{z, h) = h\g{z) - zd,g{z)]/[z - h'g{z)] . (6) 

Setting J~l, one has for small z [2^ 

GoiO,z)^ Xzln[~i^izf] Goil,z)^iiy (7) 

with A = l/VSn and fj, = v = 1/3 . The z-dependences 
of Eq. ([7]) are a direct consequence of the massless Dirac 
spectrum for the nearest neighbour tight binding model 
on the honeycomb lattice, and hold with appropriate val- 
ues for A, fi and v throughout the gapless phase. 

In this way we obtain the asymptotic be- 
haviour for small h of the vacancy energy 
£{h) ^ — /iz^[2Aln(l/ft,)]^^/^ and magnetisation 

m{h) = -dh£{h) - j/[2Aln(l//i)]~i/2 

It is apparent from the form of the projected Zeeman en- 
ergy Hz that this magnetisation is entirely localized on 
sites adjacent to the vacancy. Strikingly, each of these 
sites, labelled ri,r2,r3 in Fig. [Tfc), carries a separate 
component {mx,my and m^, respectively) with relative 
magnitude proportional to the corresponding component 
of h. From Eq. ([5]) we find the defect susceptibility 
x(/i) = dhm(h) given above, which diverges in the small 
h limit. Thermal fluctuations at a temperature T small 
compared to the exchange can be treated as fermion ex- 
citations within one flux sector, giving a linear vacancy 
susceptibility that diverges as l/[Tln(l/T)]. 

So far we have discussed response in the flux free sec- 
tor. In fact, we flnd that a vacancy has a flux bound 
to it in the ground state. Deep in a gapped phase this 
can be demonstrated analytically. Consider Jx, Jy ^ Jz- 
For Jx=Jy=0 the ground state consists of paired spins 



on 2 bonds and has degeneracy 2 . Lifting of this de- 
generacy can be studied using Jx/ Jz and Jy/Jz as small 
parameters. One obtains an effective Hamiltonian with 
leading terms proportional to the flux through each ele- 
mentary plaquette of the lattice [loj . It is minimised by 
taking zero flux through all hexagonal plaquettes but tt 
flux through the vacancy plaquette. In the gapless phase 
this leading order calculation is insufficient. Instead we 
compute the energy of different flux sectors numerically: 
a flux binding energy of — 0.027J is demonstrated by the 
negative intercept in Fig. [21 It would be intriguing to 
flnd whether mobile vacancies also bind flux, which would 
have implications for their relative statistics. 
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FIG. 2: For systems with two vacancies, the ground state 
energy difference 2AE between the sector with fluxes on va- 
cancies and the flux free sector, as a function of inverse linear 
system size 1/L. Open symbols: systems with open boundary 
conditions. Filled symbols: systems with periodic boundary 
conditions, for which there are three classes of behaviour, de- 
termined by L mod 3. Lines are quadratic fits in 1/L. Inset: 
g(z) in the presence of a flux pair separated by d = 170. 

With this conclusion in mind, we now revisit the cal- 
culation of vacancy magnetisation. The result we have 
presented for the flux free sector hinges on the behaviour 
of g{z). To carry out a similar calculation for a vacancy 
with a flux attached, we require the elements of the Green 
function for an undiluted honeycomb lattice with tt flux 
through one hexagonal plaquette, evaluated between sites 
lying on this hexagon. We flnd numerically that the flux 
generates quasi-localised modes at low energies, which 
in turn remove the singular behaviour of g{z) at small 
z. This results in a ground state vacancy magnetisation 
m{h) ~ h\n{l/h) and a low temperature linear vacancy 
susceptibility that diverges as ln(l/T). 

In detail, an efficient computation of Green function 
elements is possible for an infinite lattice with two fiuxes 
through hexagons separated by a distance d, for d < 500. 
Results of such a calculation are displayed in the inset to 
Fig. [21 where ^[(^(z)] is shown as a function of z, taken 
on the imaginary axis. It is a constant for 1 ^ |z| ^ 
1/d, but sensitive to the finite flux separation for 1/d^ 
\z\. The behaviour of the vacancy magnetisation with a 
bound flux is a consequence of constant g{z) at small z. 
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We next consider a pair of vacancies, represented by a 
tight binding model with two defect centres of the type 
shown in Fig. [He). We find that there is an interaction 
between the magnetic moments formed around each va- 
cancy, which is weak if the pair separation is large. The 
interaction sets a field scale. Above it, the susceptibility 
is a sum of two independent vacancy contributions. Be- 
low it, behaviour depends on the relative sublattices oc- 
cupied by the two defects. Moreover, since the weak-field 
response involves low-energy and long-distance features 
of the system, it is controlled by the Z2 sum of the fluxes 
associated with the vacancies, which is zero. Results for 
the flux- free sector, expressible in terms of Go, therefore 
also illustrate behaviour in the ground state sector where 
vacancies bind fluxes. 

The discussion centred on Eq. ([5]) can be repeated for 
the two-vacancy problem, with the corresponding Green 
function expressed in terms of Go, but one further matrix 
element is required in addition to the two in Eq. ([7]): the 
one between the two vacancy sites. This enters the gen- 
eralisation of £{h) and sets a scale for h. With vacancies 
on opposite sublattices and z <C |d|~^ its value is 



Go(d, 



|d|-isin(K.d-6'), 



(9) 



where K = (27r/3, — 27r/3) is the momentum at one of the 
Dirac points and 9 is the angle between d and the x axis 
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For h « h,{d) = |Go(d,0)|(ln[|l/Go(d,0)|])-i/2 
the vacancy susceptibility is large but field independent, 
while in the opposite limit (/ic(d) <^ h <^ 1) vacancies are 
independent. With both vacancies on the same sublat- 
tice, Go(d, z) vanishes as z approaches zero, and Zeeman 
energy dominates over the interaction between vacancies 
at all h. This has the remarkable consequence that a pair 
of vacancies on the same sublattice has a parametrically 
larger weak-field susceptibility {h[\n{l / h]"^/^)"^) 

than an isolated vacancy with bound flux. 

We now turn to a discussion of the influence of dis- 
order in the strength of exchange interactions. Small 
amplitude disorder in the absence of vacancies has no 
qualitative effect on some aspects of behaviour: the 
ground state has finite susceptibility and is in the flux 
free sector. Exchange randomness does however enter 
the free fermion description of states within this sector. 
We recall that without disorder this description reduces 
at low energy to two copies of the Dirac Hamiltonian, 
which in the gaplcss phase is massless Randomness 
in exchange interactions appears in the Dirac Hamilto- 
nian as a random vector potential. The consequences 
of such disorder for fermion eigenstates have been stud- 
ied extensively Most importantly in our context, 
the fermion density of states is proportional to a power 
of energy, with an exponent that depends on disorder 
strength. We introduce weak, smooth exchange disor- 
der of the form Ja{r) = J(l + ea{r)), with correlations 
(e„(r)e^(r')) = |7rA<5„^/(r - r'), where Er/W = 1 
and /(r) decreases with smoothly with r. Then for small 



A the heat capacity G has the low temperature form 
C t2/(i+^). Note that this implies that - especially 
in a system with strong spin-lattice coupling - observing 
G oc is not an experimental requirement for identify- 
ing a magnet with a gapless 'Dirac' excitation spectrum. 

In summary, we have studied the Kitaev honeycomb 
model as an example of a spin liquid that is solvable 
even in the presence of disorder. At an isolated vacancy 
a flux is bound and a local moment forms, with singular 
susceptibility in the gapless phase. Weak exchange dis- 
order does not influence the susceptibility but leads to 
a continuously variable heat capacity exponent. Finally, 
as indicated by the rich and intriguing behaviour of va- 
cancy pairs, we note that the properties of the model 
with a finite density of vacancies represents an intriguing 
and demanding open problem. 

This work was supported in part by EPSRC under 
Grant No. EP/D050952/1. 
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